Draft version February 19, 2009 

Preprint typeset using 1^1^^ style cmulatcapj v. 10/09/06 



O 

o 

(N 



6 



FORMATION, SURVIVAL, AND DESTRUCTION OF VORTICES IN ACCRETION DISKS 
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ABSTRACT 

Two dimensional hydrodynamical disks are nonlinear ly unstable to the formation of vortices. Once 
formed, these vortices essentially survive forever. What happens in three dimensions? We show with 
incompressible shearing box simulations that in 3D a vortex in a short box forms and survives just as 
in 2D. But a vortex in a tall box is unstable and is destroyed. In our simulation, the unstable vortex 
decays into a transient turbulent-like state that transports angular momentum outward at a nearly 
constant rate for hundreds of orbital times. The 3D instability that destroys vortices is a generalization 
of the 2D instability that forms them. We derive the conditions for these nonlinear instabilities to 
act by calculating the coupling between linear modes, and thereby derive the criterion for a vortex to 
survive in 3D as it does in 2D: the azimuthal extent of the vortex must he larger than the scale height of 
the accretion disk. When this criterion is violated, the vortex is unstable and decays. Because vortices 
are longer in azimuthal than in radial extent by a factor that is inversely proportional to their excess 
vorticity, a vortex with given radial extent will only survive in a 3D disk if it is sufficiently weak. This 
counterintuitive result explains why previous 3D simulations always yielded decaying vortices: their 
vortices were too strong. Weak vortices behave two-dimensionally even if their width is much less 
than their height because they are stabilized by rotation, and behave as Taylor-Proudman columns. 
We conclude that in protoplanetary disks weak vortices can trap dust and serve as the nurseries of 
planet formation. Decaying strong vortices might be responsible for the outwards transport of angular 
momentum that is required to make accretion disks accrete. 

Subject headings: accretion, accretion disks — instabilities — solar system: formation — turbulence 
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1. INTRODUCTION 

Matter accretes onto a wide variety of objects, such as 
young stars, black holes, and white dwarfs, through ac- 
cretion disks. In highly ionized disks magnetic fields are 
important, and they tri gger turbulence via the magne- 
torotational instability (jBalbus fc Hawlevlll998f) . How- 
ever, many disks, such as those arou nd young stars or 
dwarf novae, are nearly neutral (e.g.. ISano et al.l [20001 : 
iGammie fc Menoulll998D . In these disks, the fluid mo- 
tions are well described by hydrodynamics. 

Numerical simulations of hydrodynamical disks 
in two-dimensions — in the pla ne of the disk— often 
produce long-lived vor tices (iGodon fc Livio jl999t 
lUmurhan fc Regevl[200l iJohnson fc Gammie 200%. If 
vortices really exist in accretion disks, they can have 
important consequences. First and foremost, they 
might generate turbulence. Since turbulence naturally 
transports angular momentum outwards^, as is required 
for mass to fall inwards, it might be vortices that cause 
accretion disks to accrete. Second, in disks around young 
stars, long-lived vortices can trap sohd particles and 
initia te the formation of planets (iB arge fc Sommerial 
fT995h . 

Why do vortices naturally form in 2D simulations? Hy- 
drodynamical disks are stable to linear perturbations. 
However, they are nonlinearly unstable, despite some 
claims to the contrary in the astrophysical literature. 

^ CITA. Toronto, Ontario, Canada; yoram@cita.utoronto.ca 
^ Energy conservation implies that turbulence transports angu- 
lar momentum outwards; see [[3] Nonetheless, if an external energy 
source (e.g., the radiative energy from the central star) drives the 
turbulence, then angular momentum could in principle be trans- 
ported inwards. 



In two dimensions, the incompressible hydrodynamical 
equations of a disk are equiv alent to those of a non- 
rotating linear shear flow (e.g.. iLithwickl [20071 hereafter 
L07). And it has long been known that such flows 
are n onlinearly unstable (jGilll Il965t iLerner fc KnoblochI 
Il988t L07). This nonhnear instability is just a special 
case of the Kelvin-Helmholtz instability. Consider a lin- 
ear shear flow extending throughout the x-y plane with 
velocity profile v = ~qxy, where g > is the con- 
stant shear rate, so that —q is the flow's vorticity. (In 
the equivalent accretion disk, the local angular speed is 
fl — 2(7/3.) This shear flow is linearly stable to infinites- 
imal perturbations. But if the shear profile is altered by 
a small amount, the alteration can itself be unstable to 
infinitesimal perturbations. To be specific, let the alter- 
ation be confined within a band of width Ax, and let it 
have vorticity to = uj{x) (with < q), so that it in- 
duces a velocity field in excess of the linear shear with 
components Uy ~ ujAx and = 0. Then this band 
is unstable to infinitesimal nonaxisymmetric (i.e. non- 
stream-aligned) perturbations provided roughly that 



\kJ < 



q Ax 



2D instability 



(1) 



where ky is the wavenumber of the nonaxisymmetric per- 
turbation.'^ For any value of \uj\ and Ax, the band is al- 

^ More precisely, the necessary and sufficient condition for insta- 



duj I dx 



dx, where xq 



bility in the limit \uj\ < g is that l^i/l < ^ f° ^ 

is any value of x at which dw/dx = l|Gillll965l : ILerner fc KnoblochI 
119881 L07). For arbitrarily large a;, Rayleigh's inflection point the- 
orem and Fj0rtoft's theorem give necessary (though insufficient) 
criteria for instability dPrazin fc Reidi |2004I) . The former states 
that for instability, it is required that dw/dx = somewhere in 
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ways unstable to perturbations with long enough wave- 
length. Remarkably, instability even occurs when \uj\ is 
infinitesimal. Hence we may regard t his as a true non- 
linear instability. Ba lbus fc HawlevI ()2006D assert that 
detailed numerical simulations have not shown evidence 
for nonlinear instability. The reason many simulations 
fail to see it is that their boxes are not long enough in 
the y-direction to encompass a small enough non-zero 

\kyl 

In two dimensions, the outcome of this instability is 
a long-lived vortex (e.g., L07). A vortex that has been 
studied in detail is the Moore-Saffman vortex, which is 
a localized patch of spatially cons tant yorticit y super- 
imposed on a linear shear flow (S affmanlll995f ). When 
I'^l ^ 9j where to here refers to the spatially constant 
excess vorticity within the patch, and when the vortic- 
ity within the patch (w — q) is stronger than that of the 
background shear, then the patch forms a stable vortex 
that is elongated in y relative to x by the factor 



This relation applies not only to Moore-Saffman vortices, 
but also to vortices whose uj is not spatially constant. 
It may be understood as follows. A patch with char- 
acteristic excess vorticity ~ u and with Ay 3> Aa; in- 
duces a velocity field in the x-direction with amplitude 
Ux \uj\Ax, independent of the value of Ay (e.g., §6 in 
L07). As long as \uj\ < q, the y- velocity within the vor- 
tex is predominantly due to the background shear, and 
is ~ qAx. Therefore the time to cross the width of the 
vortex is tx ~ Ax/ux I/I'j-'Ii and the time to cross 
its length is ty ~ Ay /{qAx). Since these times must be 
comparable in a vortex, equation ^ follows. Equation 
^ is very similar to equation The 2D instability 
naturally forms into a 2D vortex. Futhermore, the expo- 
nential growth rate of the instability is ~ 1^1, which is 
comparable to the rate at which fluid circulates around 
the vortex. 

More generally, an arbitrary axisymmetric profile of 
U!{x) tends to evolve into a distinctive banded structure. 
Roughly speaking, bands where uj < contain vortices, 
and these are interspersed with bands where a; > 0, 
which contain no vortices. (Recall that we take the back- 
ground vorticity to be negative; otherwise, the converse 
would hold.) The reason for this is that only regions that 
have UJ < can be unstable, as may be inferred either 
from the integral criterion for instability given in foot- 
note [31 or from Fj0rtoft's theorem. For more d etail on 
vortex dynamics in shear flows, see the review bv lMarcui 
(fT99l . 

What happens in three dimensions? To date, numer- 
ical simulations of vort ices in 3D disks ha v e bee n re- 
ported in two papers. iBarranco fc Marcui (|2005f) ini- 
tialized their simulation with a Moore-Saffman vortex, 
and solved the anelastic equations in a stratifled disk. 
They found that this vortex decayed. As it decayed, new 
vortices were formed in the disk's atmosphere, two scale 
heights above the midplane. The ne w vortices survived 
for the duration of the simulation. IShen et all ()2006f ) 

the flow, i.e. that the velocity field must have an inflection point. 
ILovelace et al.l II1999I ) generalize Rayleigh's inflection point theo- 
rem to compressible and nonhomentropic disks. 



performed both 2D and 3D simulations of the compress- 
ible hydrodynamical equations in an unstratified disk, 
initialized with large random fluctuations. They found 
that whereas the 2D simulations produced long-lived vor- 
tices, in three dimensions vortices rapidly decayed. 

Intuitively, it seems clear that a vortex in a very thin 
disk will behave as it does in 2D. And from the 3D sim- 
ulations described above it may be inferred that placing 
this vortex in a very thick disk will induce its decay. 
Our main goal in this paper is to understand these two 
behaviors, and the transition between them. A crude ex- 
planation of our final result is that vortices decay when 
the 2D vortex motion couples resonantly to 3D modes, 
i.e., to modes that have vertical wavenumber ^ 0. As 
described above, a vortex with excess vorticity \uj\ has 
circulation frequency \uj\, and ky/k^ \^\/<l, where k^ 
and ky are its "typical" wavenumbers. Furthermore, it is 
well-known that the frequency of axisymmetric {ky = 0) 
inertial waves is ^k^/ ^Jk-^ + k1 (see eq. [S]). Equating 
the two frequencies, and taking the kx of the 3D mode to 
be comparable to the kx of the vortex, as well as setting 
<7 = 351/2 for a Keplerian disk, we find 

kz ^ ky (3) 
as the condition for resonance. Therefore a vortex with 
length Ay will survive in a box with height Az < Ay, 
because in such a box all 3D modes have too high a fre- 
quency to couple with the vortex, i.e., all nonzero kz ex- 
ceed the characteristic ky ~ l/Aj/. But when Az > Ay, 
there exist kz in the box that satisfy the resonance condi- 
tion fS]), leading to the vortex's destruction. This conclu- 
sion suggests that vortices live indefinitely in disks with 
scale height less than their length {h < Ay) because 
in such disks all 3D modes have too high a frequency 
for resonant coupling. This conclusion is als o consistent 
with the si mulations of IBarranco fc Marcus! (2005) and 
IShen et al] ([2006). Both of these works initialized their 
simulations with strong excess vorticity \uj\ ~ q, corre- 
sponding to nearly circular vortices. Both had vertical 
domains that were comparable to the vortices' width. 
Therefore both saw that their vortices decayed. Had 
they initialized their simulations with smaller |u;|, and 
increased the box length Ly to encompass the resulting 
elongated vortices, both would have found long-lived 3D 
vortices. IBarranco fc Marcusf s discovery of long-lived 
vortices in the disk's atmosphere is simple to understand 
because the local scale height is reduced in inverse pro- 
portion to the height above the midplane. Therefore 
higher up in the atmosphere the dynamics becomes more 
two-dimensional, and a given vortex is better able to sur- 
vive the higher it is.^ 

1.1. Organization of the Paper 

In S2]we introduce the equations of motion, and in fJ3] 
we present two pseudospectral simulations. One illus- 
trates the formation and survival of a vortex in a short 
box, and the other illustrates the destruction of a vortex 
in a tall box. 

In §^2115] we develop a theory explaining this behav- 
ior. The reader who is satisfied by the qualitative de- 
scription leading to equation ([3]) may skip those two 

* However, IBarranco fc Marcus! II2005I) also include buoyancy 
forces in their simulations, which we ignore here. How buoyancy 
affects the stability of vortices is a topic for future work. 



3 



sections. The theory that we develop is indirectly re- 
lated to the transient amplification scenario for the 
generation of turbulence. Even though hydrodynam- 
ical disks are linearly stable, linear perturbations can 
be transiently amplified before they decay, often by a 
large factor. It has been proposed that sufficiently 
amplified modes might couple nonlinear ly, leading to 
turbulence fe.g.. IChagefishvifi et all l2003l : lYeckol |2004| : 
lAfshordi et all 120051) . However, to make this proposal 
more concrete, one must work out how modes couple 
nonlinearly. In L07, we did that in two dimensions. We 
showed that the 2D nonlinear instability of equation H]) 
is a consequence of the coupling of an axisymmetric mode 
with a transiently amplified mode, which may be called 
a "swinging mode" because its phasefronts are swung 
around by the background shear. In ^we show that the 
3D instability responsible for the destruction of vortices 
is a generalization of this 2D instability. It may be under- 
stood by examining the coupling of a 3D swinging mode 
with an axisymmetric one. 3D modes become increas- 
ingly unstable as \kz\ decreases, and in the limit that 
kz ^ 0, the 3D instability matches smoothly onto the 
2D one. Thicker disks are more prone to 3D instability 
because they encompass smaller \kz\. 

2. EQUATIONS OF MOTION 

We solve the "shearing box" equations, which approx- 
imate the dynamics in an accretion disk on lengthscales 
much smaller than the distance to the disk's center. We 
assume incompressibility, which is a good approxima- 
tion when relative motions are subsonic. We also neglect 
vertical gravity, and hence stratification and buoyancy, 
which is an oversimplification. To fully understand vor- 
tices in astrophysical disks, one must consider the effects 
of vertical gravity as well as of shear and rotation. In this 
paper, we consider only two pieces of this puzzle — shear 
and rotation. Adding the third piece — vertical gravity — 
is a topic that we leave for future investigations. See also 
the Conclusions for some speculations. 

An unperturbed Keplerian disk has angular velocity 
profile 0(r) oc r~^/^. In a reference frame rotating at 
constant angular speed fio = ^(''0)7 where tq is a fidu- 
cial radius, the incompressible shearing box equations of 
motion read 

dtv + v ■Vv = -2noZXv + 2qnoxx-VP/p , (4) 

V-v^O (5) 

adopting Cartesian coordinates x, y, z, which are related 
to the disk's cylindrical r, 9 via x = r — tq and y = 
ro{9 — flot), X and z are unit vectors, and 



dlnr 



2 " 



(6) 



We retain q and fio as independent parameters because 
they parameterize different effects: shear and rotation, 
respectively. The first term on the right-hand side of 
equation ([4]) is the Coriolis force, and the second is what 
remains after adding centrifugal and gravitational forces. 
Decomposing the velocity into 

V = —qxy + It , (7) 
where the first term is the shear flow of the unperturbed 
disk, yields 

{dt — qxdy) u + u • Vm = qu^y — 2^lz Xu — VP/p(8) 
V-M==0, (9) 



dropping the subscript from ilg, as we shall do in the 
remainder of this paper. An unperturbed disk has m = 0. 

In addition to the above "velocity-pressure" formula- 
tion, an alternative "velocity-vorticity" formulation will 
prove convenient. It is given by the curl of equation ([8]), 

{dt-qxdy)u} = -qyw^, + {2n-q)d,u+Vx[u X (10) 

where 

u; = V X M (11) 

is the vorticity of u. Equation (fTO| . together with the 
inverse of equation (lll|) 

M^-V^^Vxa;, (12) 

form a complete set. 

Equation (I10|) implies that the total vorticity field is 
frozen into the fluid, because it is equivalent to 



where 



StWtot = Vx(t; X LOtot) , 
Wtot = (2ri - q)z + ijj 



(13) 
(14) 



is the total vorticity; note that ~qz is the vorticity of 
the unperturbed shear flow in the rotating frame, and 
hence (2i7 — q)z is the unperturbed vorticity in the non- 
rotating frame. The vorticity-velocity picture is similar 
to MHD, where it is the magnetic field that is frozen- in 
because it satisfies equation in place of uJtot ■ How- 
ever, in MHD the velocity field has its own dynamical 
equation, whereas in incompressible hydrodynamics it is 
determined directly from the vorticity field via equation 

m- 

3. TWO PSEUDOSPECTRAL SIMULATIONS 

The pseudospectral code is described in detail in the 
Appendix of L07. It solves the velocity-pressure equa- 
tions of motion with an explicit viscous term 

j/V^w (15) 

added to equation ([5]). (In L07, we did not include this 
term because we only considered inviscid flows.) The 
equations are solved in Fourier-space by decomposing 
fields into spatial Fourier modes whose wavevectors are 
advected by the background flow —qxy. As a result, 
the boundary conditions are periodic in the y and z 
dimension, and "shearing per iodic" in x. M ost of our 
techniques are standard fe.g..|Maron fc Goldrc ich 2001] 
lRogallol[T98lt iBarranco fc Marcusll2006f) . One exception 
is our method for remapping highly trailing wavevectors 
into highly leading ones, which is both simpler and more 
accurate than the usual method. In addition, our remap- 
ping does not introduce power into leading modes, be- 
cause a mode's amplitude has always been set to zero 
before the remap. The code was extensively tested on 2D 
flows in L07. A number of rather stringent 3D tests are 
performed in this paper. We shall show that the code cor- 
rectly reproduces the linear evolution of 3D modes (SjJ]), 
as well as the nonlinear coupling between them (fj5]). We 
also show in the present section that it tracks the various 
contributions to the energy budget, and that the sum of 
the contributions vanishes to high accuracy. 

Figures [T][4] show results from two pseudospectral sim- 
ulations. One simulation illustrates the formation and 
survival of a vortex, and the other illustrates vortex de- 
struction. In the first (the "short box"), the number of 
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Fig. 1. — Vortex Formation and Survival in a Short Box: Color depicts lUz- The initial state is unstable to vertically symmetric 
{kz = 0) perturbations, and forms into a vortex. But it is stable to 3D {k^ ^ 0) perturbations, and the evolution remains two-dimensional. 
The bottom panels show horizontal slices through the boxes in the upper panels, midway through the boxes. At time=150, the vortex 
has already formed. Only fluid with lo^ < —0.08 is shown in the middle panels to highlight the vortex, and to illustrate that surfaces of 
constant lOz remain purely vertical. At time=500, the vortex still survives. Its amplitude is slowly decaying by viscosity, which acts on 
timescale=1130. We set H = 1 and q = 3/2. The number of modes in the simulation is n^, X X = 64 X 64 X 32, and the size of the 
simulation box is {Lx , Ly , L^) = (^,1, |). In this figure, Lz is to scale relative to Ly, but Lx has been expanded by a factor of 5 for 
clarity. 




Fig. 2. — Vortex Destruction in a Tall Box: The setup is identical to the short-box simulation of Figure [T] except that the height 
Lz has been increased by a factor of 4, so that it now exceeds Ly. The resulting evolution is dramatically different. The initial state is now 
unstable not only to 2D perturbations, but to 3D ones as well. In the middle panels, surfaces of constant u!z are warped, and the evolution 
is no longer vertically symmetric. In the right panels, the flow looks turbulent. 
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Fig. 3. — Energy in the Short Box: The three contributions 
to the ene rgy budg et, E^2, AiJgijeari a^nd ABvisci a^r^ defined in 
equations l|20|l - l|22| l. E^2 initially decays, and then rises to a peak 
near t ~ 200 as nonaxisymmetric perturbations turn the axisym- 
metric mode into a vortex. Subsequently, the vortex decays due 
to viscosity. The spikiness of the evolution is due to the boundary 
conditions, as explained in the text. Also shown in the bo ttom 
panel is the error due to numerical effects, Ai?crror (eq. 1241 '). It 
is unlabelled because it is mostly obscured by ABghcai-. But it is 
nearly equal to zero everywhere, showing that the code accurately 
tracks the components of the energy budget. 

Fourier modes used is x Uy x n 



64 X 64 X 32, and 
the simulation box has dimensions — j^, Ly — 1, 
and = i . In the second (the "taU box" ) , the setup is 
identical, except that it has = 2 instead of 1/2. Both 
simulations are initialized by setting 



-0.1 cos I —X 



(16) 



In addition, small perturbations are added to long- 
wavelength modes. Specifically, labelling the wavevec- 
tors as 



{kx, ky, kz 



2tt 



Jx_ 



La, 



(17) 



with integers {jx, jy, jz), we select all modes that satisfy 
\jx\ < 3, \jy\ < 3, and \jz\ < 3, and set the Fourier 
amplitude of their ujz to 10~''e*"^, where </> is a random 
phase. But we exclude the {jx,jy,jz) = (0,0,0) mode, as 
well as {jx,jy,jz) = (±lj 0, 0), which is given by equation 
(fT6)) . Finally, we set = 1, g = 3/2, ly = IQ-'^ , and 
integration timestep dt = 1/30. 

With our chosen initial conditions, the mode given by 
equation Ijl6p is nonlinearly unstable to vertically sym- 
metric {kz — 0) perturbations, and hence it tends to 
wrap up into a vortex. From the approximate criterion 
for instability (eq. [1] ) , we see that to illustrate the wrap- 
ping up into a vortex of a mode with a small amplitude. 



Fig. 4. — Energy in the Tall Box: The initial evolution is 
almost the same as that seen in the short box (Fig. |3j. But 3D 
perturbations are unstable and force the destruction of the vortex. 
In the time interval 300 600, while the initial axisymmetric 

disturbance decays in a turbulent-like state, the value of E^2 is 
significantly larger than its initial value, and ABsijcar rises nearly 
linearly in time, corresponding to nearly constant outwards trans- 
port of angular momentum in a disk. The contribution of numeri- 
cal errors to the time-integrated energy budget, ABcrror (eq. 1241 ). 
remains small throughout. 

one must make the simulation box elongated in the y- 
direction relative to the x-scale of the mode in equation 
CH). 

In the short box (Fig. [T]), the evolution proceeds just 
as it would in two dimensions. The initial mode indeed 
wraps up into a vortex, and the evolution remains ver- 
tically symmetric throughout. Once formed, the vor- 
tex can live for ever in the absence of viscosity. But 
in our simulation, there is a slow viscous decay. The 
timescale for viscous decay across the width of the vor- 
tex is ^ ^/i^kx = 1130, taking kx = 2tt/Lx. 

In the tall box (Fig. [J), the evolution is dramatically 
different. In this case, the initial state is unstable not 
only to 2D perturbations, but to 3D (fc^ ^ 0) ones as 
well. In the middle panel of that figure, we see that 
instead of forming a vertically symmetric vortex as in 
the short box, surfaces of constant uJz are warped. By 
the third panel, the flow looks turbulent. 

Figures [3]|3] shows the evolution of the energy in these 
simulations. Projecting u onto the Navier-Stokes equa- 
tion (eq. [5] with viscosity included), and spatially aver- 
aging, we arrive at the energy equation 



dt 2 



q{UxUy) 



iy(u-\7~u) 



(18) 



after applying the shearing-box boundary conditions, 
where angled brackets denote a spatial average. The time 
integral of this equation is 

E^2 — _E„2 It^o = Ai?shoar + Ai?visc , (19) 
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where 




AE^isc = i^ (u-\/'^u)dt' . (22) 
Jo 

The pseudospectral code records each of these terms, and 
Figure [3] shows the result in the short box simulation. 
At very early times, £'„2 decays from its initial value 
due to viscosity. At the same time, the small vertically 
symmetric perturbations are growing exponentially, and 
they start to give order unity perturbations by i ~ 150, 
by which time a vortex has been formed (Fig. [1]). As 
time evolves, £'„2 gradually decays due to viscosity on 
the viscous timescale = 1130. The evolution is very spiky. 
We defer a discussion of this spikiness to the end of this 
section. 

Figure [4] shows the result in the tall box. The early 
evolution of i?„2 is similar to that seen in the short box. 
Both start with the same £'„2 , and an initial period of vis- 
cous decay is interrupted by exponentially growing per- 
turbations. But in the tall box, not only are vertically 
symmetric modes growing, but modes with ^ are 
growing as well. By t ~ 150, there is a distorted vor- 
tex that subsequently decays into a turbulent-like state. 
The energy £^^12 rises to a value significantly larger than 
its initial one, and it continues to rise until t ~ 600, 
when it starts to decay. Throughout the time interval 
300 ^ t < 600, Ai?shear rises nearly linearly in time, 
showing that {uxUy) is positive and nearly constant. 

It is intriguing that {uxUy) is positive for hundreds of 
orbits, because it suggests that decaying vortices might 
transport angular momentum outwards in disks and 
hence drive accretion. Understanding the level of the 
turbulence, its lifetime, and its nature are topics for fu- 
ture work. Here we merely address the sign of {uxUy). 
The quantity {uxUy) is the flux of y- momentum in the 
-l-x-direction (per unit mass and spatially averaged). It 
corresponds to the flux of angular momentum in a disk. 
A positive {uxUy) implies an outwards flux of angular 
momentum, as is required to drive matter inwards in an 
accretion disk. (Even though the shearing box cannot 
distinguish inwards from outwards, the sign of the an- 
gular momentum within a box depends on which side 
of the shearing box one calls inwards. Therefore, out- 
wards transport of (positive) angular momentum is well- 
defined in a shearing box.) In the shearing box, any force 
that tends to diminish the background shear flow —qxy 
necessarily transports j/-momentum in the -f x-direction. 
Hence the fact that (uxUy) > in Figure [4] shows that 
the turbulence exerts forces that resist the background 
shear, as one might expect on physical grounds. One can 
also understand why (uxUy) > from energy consider- 
ations. Since A_Evisc < 0, as may be seen explicitly by 
an integration by parts, i.e. (m-V^m) = — J2i 
equation (fT9|) may be rearranged to read 

Ai?shcar = |Ai?visc| + -E„2 — Eu2\t=o ■ (23) 

If the turbulence reaches a steady state — as it approxi- 
mately does in Figure 3] during the time interval 300 < 
t < 600 — then the last two terms in the above equation 



are nearly constant, whereas |A£'visc| increases linearly 
with time. Hence A_Eshcar must also increase. The fact 
that energy dissipation implies outwards transport of an- 
gular r nomentum is a general proper ty of accretion disks 
(e.g., iLvnden-Bell fc Pringld [l974f ). Since turbulence 
always dissipates energy, it must also transport angular 
momentum outwards. However, this argument can be vi- 
olated if an external energy source drives the turbulence, 
in which case one would have to add this energy to the 
left-hand side of equation For example, the sim- 

ulations of lStone fc Balbus. (il9.96) show that convective 
disks can transport angtdar momentum inwards when an 
externally imposed heat source drives the convection. 

Also shown in the bottom panels of Figures [3][3] is the 
integrated energy error 

Ai?crror = AiSghcar + A_Evisc + -^112 |f=o ~ (24) 

due to numerical effects, which is seen to be small. (In 
Figure [3l AiJcnor is not labelled because the curve is 
mostly obscured by AEghcav] it can be seen near t ~ 200, 
and is everywhere very nearly equal to zero.) The 
fact that AE'ciTor nearly vanishes throughout the sim- 
ulations is not guaranteed by the pseudospectral algo- 
rithm. Rather, we have chosen v to be large enough 
that the algorithm introduces negligible error into the 
energy budget. To be more precise, at each timestep in 
the pseudospectral code, modes that have \jx\ > rixji 
or \jy\ > ny/3 or \j^\ > 71^/3, where jx,y,z are de- 
fined via equation (jl7[) , have their amplitudes set to zero 
("dealiased"). This introduces an error that is analogous 
to grid error in grid-based codes. By choosing v to be 
sufficiently large, it is the explicit viscosity that forces 
modes with large k to have small amplitudes, in which 
case the dealiasing procedure has little effect on the dy- 
namics. Increasing the resolution Ux x riy x would 
allow a smaller v to be chosen — implying a larger effec- 
tive Reynolds number — while keeping the energy error 
small. 

The curves of show sharp narrow spikes every time 
interval At = 10, with width ~ 1. S imilar spikes have 
been seen in othe r simulations (jUmurha n & Regev 200||; 
IShen et al.ll2006[) . but they are stronger and narrower in 
our simulations because our simulation box is elongated. 
These spikes are due to the shearing-periodic boundary 
conditions. It is perhaps simplest to understand them by 
following the evolution in fc-space, as we shall do in ^ 
(see also LOT). But for now, we explain their origin in 
real-space. By the nature of shearing-periodic boundary 
conditions, associated with the simulation box centered 
at x = are "imaginary boxes" centered at x = jLx 
with integer j — ±1,±2,---. These imaginary boxes 
completely tile the x — y plane, and each contains a vir- 
tual copy of the conditions inside the simulation box. 
The boxes move relative to the simulation box in the 
y-direction, with the speed of the mean shear at the cen- 
ter of each box, —qjLx- Therefore, in the time interval 
At — Ly/{qLx) — 10, all the boxes line up. When this 
happens, the velocity field u that is induced by the vor- 
ticity within all the boxes (via eq. [12]) becomes large, 
because all the boxes reinforce each other, and there- 
fore £'„2 exhibits a spike. Even though the shearing- 
periodic boundary conditions that we use are somewhat 
artificial, we are confident that using more realistic open 
boundary conditions would not affect the main results 
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of this paper — and particularly not the stability of ax- 
isymmetric modes to 3D perturbations. In L07, where 
we considered 2D dynamics, we investigated both open 
and shearing-periodic boundary conditions, and showed 
explicitly that both give similar results. We also feel that 
the boundary conditions likely do not affect the level and 
persistence of the "turbulence" seen in Figure [H How- 
ever, this is less certain. Future investigations should 
more carefully address the role of boundary conditions. 

4. LINEAR EVOLUTION 

In the remainder of this paper, we develop a theory ex- 
plaining the stability of vortices seen in the above numer- 
ical simulations. We first consider the linear evolution of 
individual modes, and then proceed to show how non- 
linear coupling between linear modes can explain vortex 
stability. 

The linear evolution has been considered previ- 
ouslv (lAfshordi et all [20051 : iJohnson fc Gammid 120051 : 
iBalbus fc HawlevI 12006 ). Only two aspects of our treat- 
ment are new. First, we give the solution in terms of vari- 
ables that allow the simple reconstruction of the full vec- 
tors u) and u. And second, we give the analytic expres- 
sion for matching a leading mode onto a trailing mode 
that is valid for all ky and fc^. 

The linearized equation of motion is (cq. [lOj ) 

{dt — qxdy)(jj = —qyojx + (2f2 — q)dzU . (25) 
A single mode may be written as 

=u,(feo,t)e^[^('^0'*)]-^ , (26) 

where ko is a constant vector that denotes the wavevec- 
tor at time t — 0, and the wavevector k = k{ko,t) has 
components 

ky = koy = const (27) 
kz — koz = const (28) 
kx — kox + qtky ^ const , (29) 

so that upon insertion into equation ()25|) . the time- 
derivative of the exponential cancels the term —q xdy OJ. 
The velocity field induced by such a mode is (eq. [12] ) 



where 



u{x,t)=u{ko,t)e'^\-i^°'')y'' 



.fc X a; 



(30) 



(31) 



Figure [5] sketches the evolution of wavevectors. Ax- 
isymmetric modes {ky = 0) do not move in fc-space, as 
depicted by the sphere at (1,0,0) in Figure [5l "Swing- 
ing modes" have ky ^ 0, and their kx is time-dependent. 
Their fronts of constant phase are advected by the back- 
ground shear. Swinging modes with kx/ky < 0, as de- 
picted by the two spheres near (1,-1,0) and (1,-1,1) 
in Figure [5l have phasefronts tilted into the background 
shear, i.e., they are leading modes. As time evolves, the 
shear first swings their kx through kx = 0, at which point 
their phasefronts are radially symmetric. Subsequently, 
they become trailing modes {kx/ky > 0), and approach 
alignment with the azimuthal direction {kx/ky — > oo). 

We turn now to the evolution of the Fourier ampli- 
tudes. In the remainder of this paper we drop the hats 

Gj uj , u u . (32) 



To distinguish real-space fields, we shall explicitly write 
their spatial dependence, e.g. U3{x). 

Because a; (a;) is divergenceless, ui only has two degrees 
of freedom, which we select to be lOx and 



X- {k y. lS) 



■Jyz 



where 



k^t 



iOz 



-1- h'^ 

Ky -t- 



(34) 



Our variable u;^z is pro portional to the variable U of 
IBalbus &: HawlevI ()2006D . Adopting ujx as the second 
degree of freedom enables the full vectors to be recon- 
structed as 

kx{k X x) k X X 
" -~^yz—r (35) 



Sz 



kx{k X x) 



tv n-7/ 



lUJx 



i^yz 

k X X 



(36) 



The linearized equation (|25p is expressed in terms of 
these degrees of freedom as 



kyz d f UJ^ 
qky dt \ ^y: 



2 1+t2 





UJy 



after introducing the epicyclic frequency, 

K= ^/2Q{2n-q) , 
with K = in a Kcplcrian disks, and 

kx 



(37) 
(38) 

(39) 
(40) 



q ky 

As long as ky ^ 0, t varies in time through its depen- 
dence on kx = kox + qtky. 
For axisymmctric modes {ky = 0), t is constant and 



^UJyZ = , 



the solution of which is sinusoidal with frequency 
nkz/ \/k^~-\~k^. Axisymmctric modes with phasefronts 
aligned with the plane of the disk {kx = ky = 0) have in- 
plane fluid velocities, and they oscillate at the epicyclic 
frequency of a free test particle, k. But axisymmctric 
modes with tilted phasefronts have slower frequencies, 
because fluid pressure causes deviations from free epicy- 
cles. In the limit of vertical axisymmetric phasefronts 



{kz 



0), the effects of rotation disappear en- 



tirely, and this zero-frequency mode merely alters the 
mean shear flow's velocity profile. 

For swinging modes {ky 7^ 0), it is convenient to em- 
ploy T as the time variable. Since 



kyz d, 



qky dt dr ' 



we have 



= 



(42) 



(43) 



dr^ ' 1 + r^ 

(jBalbus &: HawlevI l2006l ). Figure [Bj plots numerical so- 
lutions of this equation, and shows that it matches the 



Fig. 5. — Evolution of Wavevectors: Modes have constant ky and k^, and kx = gtfcy+const. The three spheres depict modes that 
play important roles in nonlinear instability. The mode at (1, 0, 0) does not move in fc-space. The other two modes are swinging modes 
that are depicted in the leading phase of their swing. They will become trailing after crossing through the radially symmetric plane. The 
mode crossing through (1, —1, 0) is responsible for 2D instability that forms vortices. The one crossing through (1, —1, 1) is responsible for 
3D instability that destroys vortices. 



output from the pseudospectral code, as well as the an- 
alytic theory described below. Given LUyz, it is trivial to 
construct w and u from 



UJx 



K diVyz 



(44) 



and equations ([35)) and (1361) . 

For highly leading or trailing modes (|t| 3> 1), equa- 
tion (j43p has simple power-law solutions. 



ujyz^uJA\T\— +lub\t\— , |t| > 1 (45) 

(|Balbus &: Hawlevll2006[ ). where uja and lob are constants 
that we shall call the "normal- mode" amplitudes, and 



S= y/l - 4/32 



(46) 



which is imaginary for 1/3 1 > 1/2. As a mode's wavevector 
evolves along a line in fc-space, its amplitude is oscillatory 
if this line is much closer to the fcz axis than to the ky one, 
and non-oscillatory if the converse is true. The transition 
occurs at = 1/2. This behavior may be understood 
as a competition between shear and epicyclic oscillations. 
The timescale for to change by an order-unity factor 
due to the shear is ishcar \kx/kx\ = \kx/Qky\, and the 



timescale for epicyclic oscillations of axisymmetric modes 



is t, 



cpi 



\kx/kz\ for \kx\ » \kz\. Therefore 



ishoar/^cpi, and when |/3| ^ 1 the epicyclic time is shorter 
and so the mode's amplitude oscillates as its wavevector 
is slowly advected by the shear. But when |/3| <C 1 the 
shear changes the wavevector faster than the amplitude 
can oscillate. 

The solution (|45p breaks down in mid-swing. As 
a swinging wave changes from leading to trailing, 
its "normal-mode amplitudes" change abruptly on the 
timescale that r changes from ±1 to =Fl via 



LOA 
LJb 



Taa Tab 
Tba Tbb 



UA 

LJb 



(47) 



trail \ / \ / lead 

where the transition matrix has components 

Taa^-Tbb^ CSC [St: /2) (48) 



BA - 



cot(^7r/2)^ 



■)S+1 



11-5 r(i + V2)^ 



^9) 



Tab Sl + ST{l/2 + S/2y 

and determinant = — 1, and hence is its own inverse. The 
components are complex when |/3| > 1/2. To derive these 
components, we took advantage of t he fact that equation 
((43)) has hypergeometric solutions ()Johnson fc Gammid 
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Fig. 6. — Linear Evolution of Mode Amplitudes for Three 
Values of /3: Time runs from right to left. Solid curv es show the 
exact, numerically integrated solution of equation (|43j. The initial 
value of dwyz/dr was chosen so that iog = initially (eq. 1451 ). 
Dashed lines show the analytic solution (eq. [45]) with constant uja 
and ujg = for r > 0; while for r < 0, the normal mode a mpli tudes 
are set to different constants that are given by equation (1471 . We 
exclude the domain |t| < 1 from the dashed curve, because the 
analytic approximation does not apply there. Circles show output 
from the pseudospectral code, integrated with a timestep dt = 1/15 
and with the viscosity set to zero. 



120051 : iBalbus fc Hawlevll2006f ). and matched these onto 
the normal-mode solution given above. We omit the un- 
enlightening details. 

5. NONLINEAR EVOLUTION: FORMATION AND 
DESTRUCTION OF VORTICES 

5.1. Qualitative Description 

The instability that destroys vortices is a generaliza- 
tion of the one that forms them. We review here how 
vortices form, before describing the instability that de- 
stroys them. In §5.2( we make this description quantita- 
tive. 

Vortices form out of a nonlinear instability that in- 
volves vertically symmetric [kz — 0) modes. (See L07 
for more details of the 2D dynamics than are presented 
here.) Consider the two vertically symmetric modes 
shown in Figure O the "mother" mode at (1,0,0) and 
the "father" mode that is depicted crossing through 
(1,-1,0). Triplets of integers {jx,jy,jz) label values 
of wavevectors (kx, ky, kz) (for example, via [17] ). The 
mother is both axi- and vertically-symmetric, and the 
father is a leading swinging mode. 

As the father swings through radial symmetry, i.e. 
as it crosses through the point (0,-1,0), its veloc- 
ity field is strongly amplified by the background shear. 
This can be seen from fj4l which shows that swing- 
have ujyz=C0Tist., and hence 



ing modes with k^ 

Ux = i[U}yz/kyz)/{l 



T crosses through 0. When the father is near the peak 
of its transient amplification (|r| < 1), it couples most 
strongly with the mother, and they produce a "son" near 
(1, -1, 0) = (1, 0, 0) + (0, -1, 0). The son will then swing 
through radial symmetry where it will couple (oedipally) 
with the mother to produce a grandson near (1,-1,0), 
which can repeat the cycle. We summarize this 2D in- 
stability feedback loop as 

linear amplification : (1,-1,0) (0,-1,0) 
nonlinear coupling : (0,-1,0) + (1,0,0) (1,-1,0) 

The criterion for instability is simply that the amplitude 
of the son's ojyz be larger than that of the father. As 
shown in L07, if instability is triggered, its nonlinear out- 
come in two dimensions is a long-lived vortex. 

The three-dimensional instability that is responsible 
for destroying vortices is a straightforward generaliza- 
tion. The mother mode is still at (1,0,0), but now the 
father mode starts near (1,-1,1). Symbolically, the feed- 
back loop is 



linear amplification : 
nonlinear coupling : (0, 



(1,-1,1)^(0,-1,1) 
-1,1) + (1, 0,0)^(1, -1,1) 



The 2D instability described above is just a special case 
of this 3D one in the limit that kz = 0. In general, the 
stability of a mother mode at (1, 0, 0) with given k^ = kx 
and ijj = u) depends on the ky and kz of the father- 
mode perturbations (as well as on the parameters q and 
SI). Which ky and kz are accessible in turn depends 
on the dimensions Ly x Lz of the simulation box — or 
equivalently, on the circumferential distance around a 
disk and the scale-height. In ij5.2i we map out quan- 
titatively the region in the ky — k^ plane that leads to 
instability. For now, it suffices to note that the unstable 
region has jfej;! < \kxCj\/q and Ifc^j < \ky\. We conclude 
that a given mother mode suffers one of three possible 
fates, depending on Ly and Lz- 



1. li Ly is less than a critical value (~ ql\ujkx\), then 
the mother mode is stable to all perturbations. 



T^), which becomes largest when 



2. If Ly is larger than this critical value, then the 
mother mode is unstable to vertically symmetric 
{kz — 0) perturbations; if in addition Lz is suffi- 
ciently small that all modes with kz ^0 are stable, 
then the mother mode turns into a long-lived vor- 
tex (Figure [T]). 

3. If both Ly and L^ are sufficiently large, the mother 
mode is unstable both to vertically symmetric and 
to 3D perturbations. When this happens, the 
mother starts to form a vortex, but this vortex is 
3D-unstable. The result is turbulence (Figured]). 

There is also a possibility that is intermediate between 
numbers 2 and 3: if the conditions described in number 
2 hold, the essentially 2D dynamics that results can non- 
linearly produce new mother modes that are unstable 
to 3D perturbations. In this paper, we shall not con- 
sider this possibility further, since it did not occur in the 
pseudospectral simulations of f}3j We merely note that 
in our simulations of this possibility (not presented in 
this paper), we found that when the new mother modes 
decayed, they also destroyed the original mother mode. 
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Fig. 7. — Nonlinear Evolution of 3D Instability: Time runs 
from right to left. In t he tw o le ft panels, lines show numerical solu- 
tions of equations l|50| l and ((52}, as well as the grandson's equation. 
Also shown as circles are the output from a pseudospectral simula- 
tion, showing excellent agreement with the "exact" solutions. The 
following parameters have been chosen: H = 1, g = 3/2, Cb = 0.005, 



-27r/30, fcz = 0.45|fca| =^ /9 = -0.3. The small 



disagreement between pseudospectral and exact solutions for ui'^ 
at r < — 2f is due to the conju gate mod es th at, for simplicity, we 
have not included in equations I I50I I and I I52I I ; see footnote |5] The 
two right panels show the mode amplitudes, defined via equation 
IIAll l for the son, and similarly for the father and grandson. Al- 
though these two right panels contain the same information as the 
left ones, they are helpful in constructing the analytic form of the 
growth factor x (see App endix). With the parameters chosen for 
this figure, equation 1 154 I I predicts x = —2.2 for the amplification 
factor between successive generations, in agreement with that seen 
in the figure. 

5.2. The Stability Criterion 

To quantify the previous discussion, we choose an ini- 
tial state as in Figure[5l with the mother mode at (1, 0, 0) 
and the father a leading mode crossing through (1,-1,1). 
The son mode, not depicted in the figure, is initially 
crossing through the point (2,-1,1). We set its ini- 
tial vorticity — as well as the initial vorticity of all modes 
other than the mother and father — to zero.^ The father's 
wavevector and Fourier amplitude are labelled as in SjH 

^ We ignore the complex conjugate modes for simplicity. Since 
aj(£c) is real-valued, each mode with wavevector and amplitude 
(fc,aj) is accompanied by a conjugate mode that has (— fc, a;*). In 
our initial state, there are really four modes with non-zero am- 
plitudes: the mother at (1,0,0) and its conjugate at (—1,0,0), 
and the father and its conjugate. We may ignore the conjugate 
modes because they do not affect the instability described here. 
As shown in L07 for the 2D case, their main effect is that when 
the son swings through (0, —1, 1), not only does it couple with the 
mother at (1, 0, 0) to produce a grandson at (1, —1, 1), but it also 
couples with the conjugate mother at (—1, 0, 0) to partially kill its 
father, which is then at (—1, —1, 1) (bringing to mind the story of 
Oedipus). But since the father is a trailing mode at this time, it 
no longer participates in the instability. Nonetheless, the conjugate 
modes do play a role in the nonlinear outcome of the instability. 



and the mother's and son's are labelled with bars and 
primes: 

father : k , uj 

mother : k = kx , u> = cuz 

son : k' = kx + k , us' 

Note that fc=const., and u)x — ^ because the vorticity 
must be transverse to the wavevector. We also set Qy = 
0; otherwise Uz ^ 0, which corresponds to a mean flow 
out the top of the box and in through the bottom. 

At early times, the father mode swings through the 
point (0, —1, 1). Since the only other nonvanishing mode 
at this time is the mother, there are no mode cou- 
plings that can nonlinearly change the father's ampli- 
tude. Therefore its amplitude is governed by the linear 
equation (|37p . which we reproduce here as 



1 1 

2 TFT+T' 





(50) 



During its swing, it couples with the mother to change 
the amplitude of the son. The linear part of the son's 
evolution is given by the above equation with primed 
vorticity and wavevector in place of unprimed. The non- 
linear part is given by 

^u)' = ik' X {u X UJ + u X uj) (51) 

dt nonlin 

(eq. [TU]) where u — —i{uj/k)y and u — ik X uj/k^ 
(eq. [31]). Adding the linear and nonlinear parts, and 
re-expressing in terms of our chosen degrees of freedom, 
we find 



d 

d^ 




2 



1 K'' 1 

"2?FT+(7+?p' 




^x 
L0[. 



1 + 



LOx 
LU„, 



(52) 



(53) 
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where the dimensionless constant 

kx 

depends on both the mother's and father's wavevectors. 
It is the father's r = k^/ky^ that is being used as the 
time-coordinate for evolving the son's amplitude. The 
grandson's equation is the obvious extension: denoting 
the grandson's amplitudes with double primes, one need 
only make the following replacements in equation (j52p : 
u)' u)", UJ uj' , and t —t t + f. Subsequent genera- 
tions evolve analogously. 

The father's equation (I50|) is easily solved, as shown in 
21 Inserting this solution into equation ([5^ produces a 
linear inhomogeneous equation for the son's amplitude, 
and similarly for the grandson's. Figure [7] plots numeri- 
cal solutions of these equations. Also shown as circles are 
output from a pseudospectral simulation, showing excel- 
lent agreement. 

In the Appendix, we solve equation (|52p analytically 
to derive the amplification factor x, which is the ratio of 
the son's amplitude at any point in its evolution (e.g., 
when it is radially symmetric), to the father's amplitude 
at the same point in its evolution. We find 

-1 + 5 , qVL,, ,A r(l + (5/2) 



X = — T 
1 



(52 



^ + f(^-^)yr(i/2- 



5/2) 
(54) 
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Fig. 8.— Curves of Marginal Stability for a Mother Mode With Q = 0.005 (left panel) and d> = 0.05 (right panel): Left 
panel corresponds to Figurc[7|and right panel corresponds to the pseudospectral simulations of Figures lH2l We set f2 = 1 and q = 3/2. To 
make the solid lines in these plots (the "exact solutions"), we repeated the integrations that produced the lines in Figure [T] but varying 
fcz for eac h k y until perturbations neither grew nor decayed. The dashed line in the left panel shows that the analytic approximation of 
equation H54| l agrees reasonably well with the exact solution. We do not show equation H54| l in the right panel because the agreement 
is poorer there. Right panel shows two X's for the values of the smallest non-zero \ky\ and \kz.\ in the simulations of Figures IH2I i.e., 
\kylk\ = Lx/Ly = 0.067 for both simulations, and |fcz/A:| = L^/L^ = 0.13 for the short box and 
contains a 3D-unstable mode that leads to the destruction of the vortex into a turbulent-like state, 
and is stable to 3D perturbations. 



= 0.033 for the tall box. The tall box 
The short box contains no such mode, 



where 5 — -y/l — 4/3^. Equation is apphcable in 
the limit \u)\/q < 1. For 2D modes (/? = ^ (5 = 1), 
it recovers equation 42 of L07 (see also eq. [1] of this 
paper): 

uj k 

X2D=~TT-—. (55) 

q h,y 

Marginally stable modes have \x\ = 1. Figure [5] plots 
curves of marginal stability. The left panel is for the 
case UJ = 0.005, as in Figure [71 and the right panel is for 
to = 0.05, as in the pseudospectral simulations presented 
at the outset of this paper (eq. [16]; Figs. [T][2]). The left 
panel shows that equation (|54p gives a fair reproduction 
of the exact curve. We do not show equation ([54]) in 
the right panel, because it gives poorer agreement there 
(since |(x)|/g is too large). In the right panel we also 
plot X's for the values of the smallest nonvanishing 3D 
wavenumbers in the simulations of Figures [TJH In the 
short-box simulation, all 3D modes lie in the stable zone. 
Therefore the dynamics remains two-dimensional. But 
in the 3D box, there is a 3D mode in the unstable zone 
that destroys the vortex and gives rise to turbulent-like 
behavior. 

It is interesting to consider briefly how the instability 
described here connects with the Rayleigh-unstable case, 
which occurs when < . At small \ky\, the marginally 
stable curves in Figure [8] are given by = 1/2, where 
/? = {K/q){kz/ky). Hence if one decreases k from its 
Keplerian value the marginally stable curve becomes 
steeper in the k^ — ky plane, and an increasing number 
of 3D modes become unstable. As k — > 0, if a 2D mode 
with some ky is unstable, then so are all 3D modes with 
the same ky. Therefore any 2D-unstable state is also 
3D-unstable, and any forming vortex would decay into 
turbulence. 

6. CONCLUSIONS 



Our main result follows from Figure [3 which maps 
out the stability of a "mother mode" (i.e., a mode with 
wavevector kx and amplitude cj) to nonaxisymmetric 3D 
perturbations. A mother mode is unstable provided that 
the ky and k^ of the nonaxisymmetric perturbations sat- 
isfy both \ky\ < kuj/q and \kz\ ^ \ky\, dropping order- 
unity constants. Based on this result, we may understand 
the formation, survival, and destruction of vortices. Vor- 
tices form out of mother modes that are unstable to 2D 
(fcz — 0) perturbations. Mother modes that are unsta- 
ble to 2D modes but stable to 3D (fc^ ^ 0) ones, form 
into long-lived vortices. Mother modes that are unstable 
to both 2D and 3D modes are destroyed. Therefore a 
mother mode with given fc and Co will form into a vortex 
if the disk has a sufficiently large circumferential_ extent 
and a sufficiently small scale height, i.e., if r > k~^q/uj 
and h < fc^^g/tD, where r is the distance to the center 
of the disk, and h is the scale height. Alternatively, the 
mother mode will be destroyed in a turbulent-like state 
if both r and h are sufficiently large (r > k~^q/uj and 
h > k-^q/uj). 

Our result has a number of astrophysical consequences. 
In protoplanetary disks that do not contain any vor- 
tices, solid particles drift inward. Gas disks orbit at 
sub-Keplerian speeds, Ugas ~ VLr{l — rj), where ilr is the 
Keplerian speed and 77 ~ (cs/f2r)^, with Cs the sound 
speed. Since solid particles would orbit at the Keple- 
rian speed in the absence of gas, the mismatch of speeds 
between solids and gas produces a drag on the solid par- 
ticles, removing their angular momentum and causing 
them to fall into the star. For example, in the min- 
imum mass solar nebula, meter-sized particles fall in 
from 1 AU in around a hundred years. This rapid infall 
presents a serious problem for theories of planet forma- 
tion, since it is difficult to produce planets out of dust in 
under a hundred years. Vortices can solve this problem 
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(|Barge &: Sominerialll995D . A vortex that has excess vor- 
ticity —u) and radial width 1/fc can hah the infaU of par- 
ticles provided that ui/k > {^r)r], because the gas speed 
induced by such a vortex more than compensates for the 
sub-Keplerian speed induced by gas pressure.® Previ- 
ous simulations implied that 3D vortices rapidly decay, 
and so cannot prevent the rapid infall of so lid particles 
(jBarranco fc Marcusll200l IShen et al.ll2006[ ). Our resuh 
shows that vortices can survive within disks, and so re- 
stores the viability of vortices as a solution to the infall 
problem. 

A more important — and more speculative — 
application of our result is to the transport of angular 
momentum within neutral accretion disks. In our 
simulation of a vortex in a tall box, we found that as 
the vortex decayed it transported angular momentum 
outward at a nearly constant rate for hundreds of 
orbital times. If decaying vortices transport a significant 
amount of angular momentum in disks, they would 
resolve one of the most important outstanding questions 
in astrophysics today: what causes hydrodynamical 
accretion disks to accrete? To make this speculation 
more concrete, one must understand the amplitude and 
duration of the "turbulence" that results from decaying 
vortices. This is a topic for future research. 

In this paper, we considered only the effects of rotation 
and shear on the stability of vortices, while we neglected 
the effect of vertical gravity. There has been a lot of 
research in the geophysical community on the dynamics 



of fluids in the presence of vertical gravity, since stably 
stratified fluids are very common on Earth — in the atmo- 
sphere, oceans, and lakes. In numerical and laboratory 
experiments of strongly stratified flows, thin horizontal 
"pancake vortices" often form, and fully developed tur- 
bulence is characterized by thin horizontal layers, (e.g., 
iBrethouwer et al.ll2007D . Pancake vortices are stabihzed 
by vertical gravity, in contrast to the vortices studied 
in this paper which are stabilized by rotation. Gravity 
inhibits vertical motions because of buoyancy: it costs 
gravitational energy for fluid to move vertically. The re- 
sulting quasi-two-dimensional flow can form into a vor- 
tex. We may speculate that in an astrophysical disk 
vertical gravity provides an additional means to stabi- 
lize vortices, in addition to rotation. But to make this 
speculation concrete, the theory presented in this paper 
should be extended to include vertical gravity. 

We have not addressed in this paper the origin of the 
axisymmetric structure (the mother modes) that give rise 
to surviving or decaying vortices. One possibility is that 
decaying vortices can produce more axisymmetric struc- 
ture, and therefore they can lead to self-sustaining tur- 
bulence. This seems to us unlikely. We have not seen 
evidence for it in our simulations, but this could be be- 
cause of the modest resolution of our simulations. Other 
possibilities for the generation of axisymmetric structure 
include thermal instabilities, such as the baroclinic in- 
stability, or convection, or stirring by planets. This, too, 
is a topic for future research. 



° We implicitly assume here that the stopping time of the par- 
ticle due to gas drag is comparable to the orbital time, which is 
true for meter-sized particles at 1 AU in the minimum mass solar 
nebula. A more careful treatment shows that a vortex can stop 
a particle wit h stopping time ts provided that O/k > {Qr){Qts)r] 
(lYoudin"200a). 

^ Billant fc Chomad ((2003) show that a vertically uniform vortex 
column in a stratified (and non-rotating and non-shearing) fluid 
suffers an instability (the "zigzag instability") that is character- 
ized by a typical vertical lengthscale Xz ~ U/N, where U is the 
horizontal speed induced by the vortex, is the Brunt- Vaisala 



frequency, and the horizontal lengthscale of the vortex Lf^ is as- 
sumed to be much greater that Xz (hence the pancake structure). 
We may understand Billant & Chomaz's result in a crude fashion 
with an argument similar to that employed in the introduction to 
explain the destruction of rotation-stabilized vortices: since the 
frequency of buoyancy waves is Nkx/kz (when \kx\ "C \kz\), and 
since the frequency at which fluid circulates around a vortex is 
U / ~ kxU , there is a resonance between these two frequencies 
for vertical lengthscale 1/kz ~ U/N. 



APPENDIX 

ANALYTIC EXPRESSION FOR GROWTH FACTOR x 



In this Appendix, we derive equation ()54|) by analytically integrating equation ([52]) for the son's vorticity, given the 
father's vorticity as a function of time (|J4]), and taking the mother's vorticity lu to be constant, which is valid when 
the father's amplitude is small relative to the mother's. The numerical integral of equation (|52p is shown in Figure [T) 



Recall that initially t = t > and r decreases in time, and typically r 3> 1. 
the son's "normal- mode" amplitudes and uj'g, defined from lo'^ and cUy^ via 



where 

Substituting this into equation 
equation if we approximate 1 4 
term produces 

d , (D 1 , ,, 1+.S 



K 1 — 5 1 
2/30 2r' I 



K 1 + S I 

1.^'! - 



It simplifies the analysis to work with 
eqs. [44] and [45] ^ 



(Al) 




t' = T + f (A2) 
the time derivative of the above matrix cancels the homogeneous term in that 
i; r ^, which holds until just before the time that t — ~f. The inhomogeneous 



dT 



1 + 6 1 



2i3n 1 



> -i 



(A3) 



Since lu^ and ujyz are known O, a straightforward integration yields uj'^ just before r = — r. To perform this integral, 
we resort to some approximations, guided by the solution shown in Figure [7] For the first term, we need 



(r + f)- 



1 



dr^ 



(r + r) 2 dr + LOB 



'-{t + f)~dT 



(A4) 
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_xl+S -l + S _I -3 + S , . , 1 + 

1-'-' 1 , . _ ^— L' I n1 /T n\o. 

2/32 



''^bt'^^" -^sf* / s— (l-s) — (A5) 



^ ^l + <5 r(l/2 + V2) 

--^-2 r(l + ^/2) 

where e is a parameter that satisfies 1 ^ e ^ 1/f. In the first fine, we used equation ([^5]) . and we discarded the loa 
mode from the second integral because the lob mode increases faster with increasing |r|. From Figure [7l lo'^ nearly 

vanishes until r ~ 0. Therefore, in the second line we approximated the first integral as —{f^^ / 0^)du>yz / dT\-e. The 
third line holds in the limit of small e. The other three terms in equation (|A3p are integrated similarly, yielding 



,/ ,. ;- r(i/2 + V2) / qn \ {i + 6f 

just before the time when r = — f, i.e., just before the son is radially symmetric. At this time. Figure [7| shows that 
bj'g very nearly vanishes. Therefore, just after the son is radially symmetric, it will have Lo'g = Tba^'a ("^q- |47|). with 
lo'a given by equation (IA7p . This gives the amplification factor x = ^'bI'^b that is displayed in equation 
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